Instanton Project Overview

The problem

• Consider a network with a significant number of wind farms.
• Suppose the demand forecast is accurate and generation dispatch does not change (though generators do perform droop response to handle mismatches).
• Now suppose the wind forecast is not entirely accurate. Let each wind farm's active power output be an optimization variable.
• Now for the optimization:
  • The objective is to minimize the normed distance between the vector of wind optimization variables and the wind forecast; we reason that the most concerning scenarios are those closest to the forecast, because they are most likely to occur. We might also insert a "covariance" matrix into the objective function to accommodate situations where a number of wind farms are likely to experience the same wind variations.
  • The set of constraint functions is as you expect: power balance constraints, an angle reference, a constraint on $\alpha$ to make sure all slack is taken up.
  • We add one additional constraint to force a line to its active power flow limit in one direction. After optimizing with the first line constraint in place, we reverse its direction and repeat. After that optimization is complete, we move on to the next line until all have been considered.
• In the end, we have an instanton candidate for each line constraint in the network (including both directions). Each candidate has a "score" -- its distance from the wind forecast -- and we call the one with the lowest score the instanton.

Previous work

• The instanton problem with demand uncertainty
• The nonlinear AC instanton problem
• The instanton problem with proportional generator response

Mathematical Formulation

\begin{align} \min & \frac{1}{2} \left( \rho - \rho_0 \right)^\top \left( \rho - \rho_0 \right) \\ & \text{subject to:} \\ \sum_k( Y_{ik} \theta_k) & = G_{i} + \rho_i - D_{i} \quad \forall i \in N \\ G &= G_{0} + k\alpha \\ \theta_i - \theta_k & = x_{ik} P_{lim,ik} \quad \text{for each } (i,k) \in G \\ \theta_i & = 0 \quad \text{ where bus $i$ is the angle reference} \\ \sum_i k_i & = 1 \end{align}

Notation (in order of appearance)

• $\rho_i$ is renewable generation at bus $i$ in per unit; thus, $\rho$ is a vector of decision variables.
• $R_p$ is the renewable generation forecast vector.
• $Y_{ik}$ is the $(i,k)$th element of the admittance matrix $Y$, which assumes resistances throughout the network are zero.
• $\theta_k$ is the phase angle at bus $k$.
• $G_{p,i}$ is the conventional active power generation at bus $i$.
• $D_{p,i}$ is the active power demand at bus $i$.
• $N$ is the set of buses (nodes).
• $k$ is the vector of participation factors for conventional generators (the case where $k_i=1$ corresponds to generator $i$ taking up all slack).
• $\alpha$ is the participation coefficient (mismatch), defined as

$$\alpha:= \sum D - \sum \rho - \sum G_0 ~.$$

• $G_{p0}$ is the base-case conventional active power generation vector.
• $x_{ik}$ is the reactance of line $(i,k)$.
• $P_{lim,ik}$ is the active power limit of line $(i,k)$.
• $G$ is the set of edges (lines).

A simple example

Suppose we have a triangle-shaped network with three nodes and three edges.

Two of the nodes have wind farms and the third has a conventional generator.

Assume:

• DC power flow approximation holds for this network,
• the conventional generator takes up all slack, and
• that the conventional generator node is the network's angle reference.

Variables

Our instanton optimization problem has two decision variables: $R_{1}$ and $R_{2}$.

These variables influence the two independent bus phase angles, $\theta_1$ and $\theta_2$.

Power balance and constraints

The power balance equations are:

$$\begin{align} R_{1} - d_1 &= y_{12}(\theta_1 - \theta_2) + y_{1S}\theta_1 \\ R_{2} - d_2 &= y_{12}(\theta_2 - \theta_1) + y_{2S}\theta_2 \\ \end{align}$$

Now suppose each line in the system has the same power flow limit, which we will represent as $P^{lim}$. This places a constraint on the active power flow on each line:

$$\begin{align} |y_{12}(\theta_1 - \theta_2)| &\leq P^{lim} & \implies |\theta_1 - \theta_2| &\leq \frac{P^{lim}}{y_{12}} \\ \\ |y_{1S}(\theta_1)| &\leq P^{lim} & \implies |\theta_1| &\leq \frac{P^{lim}}{y_{1S}} \\ \\ |y_{2S}(\theta_2)| &\leq P^{lim} & \implies |\theta_2| &\leq \frac{P^{lim}}{y_{2S}} \end{align}$$

Plugging in numbers...

Let's plug in some numbers and plot this set of inequalities.

Plim = 1;
y12 = 0.6;
y1S = 0.3;
y2S = 0.2;

R1 = 0.3;
R2 = 1.4;

Feasible region of optimization

Solution via KKT Conditions

To write a Lagrangian, we must first express the problem in terms of matrices and vectors. Let the set of nodes $N$ be partitioned into $N_r$, the set of nodes with renewable energy, and $N_c$, the set of nodes with conventional generation and demand. Sorting the rows of the admittance matrix $Y$, we obtain $$\begin{align} Y & = \begin{bmatrix} Y_r \\ Y_c \end{bmatrix} ~. \end{align}$$ Partitioning the remaining network parameters in this way, we can re-write the power balance constraint as follows: $$\begin{align} \begin{bmatrix}Y_r \\ Y_c \end{bmatrix}\theta + \begin{bmatrix}D_r - G_r - \rho \\ D_c - G_c \end{bmatrix} & = 0~, \end{align}$$ where $D_r$ and $G_r$ contain all components of the demand and conventional generation vectors associated with nodes that have renewable energy, and $D_c$ and $G_c$ contain components associated with nodes that have no renewable energy.

To capture conventional generation response, we replace $G$ with the following expression: $$\begin{align} G & = G_0 + k\alpha ~, \end{align}$$ where $k$ is a predetermined column vector of generator participation factors.

Each generator's droop characteristic output is given by: $$\begin{align} G^i &= G_0^i +k^i\alpha \end{align}$$ If the sum of all elements of $k$ is 1, we can sum over all $i$ and obtain the following: $$\begin{align} \sum_i G^i &= G_0^{tot} +\alpha \end{align}$$ The variable $\alpha$, then, balances total generation with total load: $$\begin{align} G_0^{tot} +\alpha + \rho^{tot} - D^{tot} &= 0 \\ \implies \alpha &= D^{tot} - \rho^{tot} - G_0^{tot} \end{align}$$ The Lagrangian term associated with $\alpha$, then, is: $$\begin{align} \lambda_\alpha\left[D^{tot} - \rho^{tot} - G_0^{tot} - \alpha \right] \end{align}$$

The slack bus angle reference is established as follows: $$\begin{align} s_{ref}^\top\theta & = 0~, \end{align}$$ where $s_{ref}$ is a column of zeros with a one in the $i$th position (bus $i$ being the angle reference).

Finally, we may saturate line $(i,k)$ in the network using the following constraint: $$\begin{align} s_{ik}^\top\theta - P_{lim,ik} & = 0 ~, \end{align}$$ where $s_{ik}$ is a column vector with positive $Y_{ik}$ in the $i$th position and negative $Y_{ik}$ in the $k$th position.

Now we have a quadratic objective function with a set of homogeneous constraints, all expressed in terms of matrices and vectors. This gives rise to the Lagrangian: $$\begin{align} \mathcal{L} \left(\rho ,\theta ,\alpha,\lambda_r,\lambda_c,\lambda_\alpha,\lambda_{ref},\lambda_{lim} \right) &= \frac{1}{2}{\left(\rho -{\rho }_{0}\right)}^{\top }\Lambda \left(\rho -{\rho }_{0}\right)+ \\ & \qquad \lambda_r^\top \left[{Y}_r \theta + D_r - (G_{0,r} + k_r\alpha) - \rho\right] + \\ & \qquad \lambda_c^\top \left[ Y_c \theta + D_c - (G_{0,c} + k_c\alpha) \right] + \\ & \qquad \lambda_\alpha\left[ D^{tot} - \rho^{tot} - G_0^{tot} - \alpha \right] + \\ & \qquad \lambda_{ref}(s_{ref}^\top\theta) + \\ & \qquad \lambda_{lim}(s_{ik}^\top \theta - P_{lim,ik}) \end{align}$$

The KKT conditions tell us to set all partial derivatives of the Lagrangian equal to zero:

$$\begin{align} \frac{\partial \mathcal{L}}{\partial \rho} = 0 &= (\rho - \rho_0)^\top \Lambda - \lambda_r^\top - \lambda_\alpha \\ \frac{\partial \mathcal{L}}{\partial \theta} = 0 &= \begin{bmatrix} \lambda_r^\top & \lambda_c^\top \end{bmatrix} \begin{bmatrix} Y_r \\ Y_c \end{bmatrix} + \lambda_{ref}s_{ref}^\top + \lambda_{lim} s_{ik}^\top \\ \frac{\partial \mathcal{L}}{\partial \alpha} = 0 &= \begin{bmatrix} \lambda_r^\top & \lambda_c^\top \end{bmatrix} \begin{bmatrix} -k_r \\ -k_c \end{bmatrix} - \lambda_\alpha \\ \frac{\partial \mathcal{L}}{\partial \lambda_r} = 0 &= Y_r \theta + D_r - (G_{0,r} + k_r\alpha) - \rho \\ \frac{\partial \mathcal{L}}{\partial \lambda_c} = 0 &= Y_c\theta + D_c - (G_{0,c} + k_c\alpha) \\ \frac{\partial \mathcal{L}}{\partial \lambda_\alpha} = 0 &= D^{tot} - \rho^{tot} - G_0^{tot} - \alpha \\ \frac{\partial \mathcal{L}}{\partial \lambda_{ref}} = 0 &= s_{ref}^\top\theta \\ \frac{\partial \mathcal{L}}{\partial \lambda_{lim}} = 0 &= s_{ik}^\top \theta - P_{lim,ik} \end{align}$$

This system of equations may be expressed in matrix form as follows: $$\begin{align} \begin{bmatrix} \Lambda & 0 & 0 & -I & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & Y_r^\top & Y_c^\top & 0 & s_{ref} & s_{ik} \\ 0 & 0 & 0 & -k_r^\top & -k_c^\top & -1 & 0 & 0 \\ -I & Y_r & -k_r & 0 & 0 & 0 & 0 & 0 \\ 0 & Y_c & -k_c & 0 & 0 & 0 & 0 & 0 \\ -\mathbf{1}^\top & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & s_{ref}^\top & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & s_{ik}^\top & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \rho \\ \theta \\ \alpha \\ \lambda_r \\ \lambda_c \\ \lambda_\alpha \\ \lambda_{ref} \\ \lambda_{lim} \end{bmatrix} & = \begin{bmatrix} \Lambda \rho_0 \\ 0 \\ 0 \\ G_{0,r} - D_r \\ G_{0,c} - D_c \\ G_0^{tot} - D^{tot} \\ 0 \\ P_{lim,ik} \end{bmatrix} \end{align}$$

---

What we've covered so far

We looked at the instanton project. This work attempts to answer the research question, "What is the smallest deviation from a wind forecast that will cause congestion in a transmission network with significant wind energy penetration?" After a high-level discussion, we formulated the problem mathematically and derived its solution in terms of the KKT conditions.

Where to next?

Next we will see how the instanton problem is formulated and solved in Julia using the JuMP mathematical programming package. We will also visualize solutions using Graphviz.